Examples
Reading Filterbanks
First, let's grab some filterbank data, we can find one from Berkeley's SETI Research Center. Julia's download function will grab the file from a URL and give us the filename back.
filename = download("https://github.com/UCBerkeleySETI/breakthrough/blob/master/GBT/filterbank_tutorial/blc04_guppi_57563_69862_HIP35136_0011.gpuspec.0002.fil?raw=true")"/tmp/jl_5EO3ZV"Now, we can read it into a Filterbank. By default, it will read all the time samples, we can however supply start and stop arguments if we want to look at a limited chunk. This is usefull as these files can get quite large.
fb = Filterbank(filename)A Filterbank file with 272 time samples
----------
az_start : 0.000000
foff : -0.002861
nifs : 1
za_start : 0.000000
nbeams : 1
machine_id : 20
tsamp : 1.073742
nbits : 32
src_raj : 71550.064000
src_dej : 471420.040000
fch1 : 1501.463413
rawdatafile : guppi_57563_69862_HIP35136_0011.0000.raw
ibeam : 1
tstart : 57563.808588
data_type : 1
telescope_id : 6
source_name : HIP35136
nchans : 65536
The data is stored as a DimArray from the DimensionalData package in the data attribute. This provides zero-cost abstractions for named indexing and fast index lookups. Our filterbank data has axes Freq and Ti, which we can use for indexing.
fb.data272×65536 DimArray{Float32,2} with dimensions:
Ti Sampled 0.0:1.0737418239999998:290.984034304 ForwardOrdered Regular Points,
Freq Sampled 1501.4634132385254:-0.00286102294921875:1313.9662742614746 ReverseOrdered Regular Points
3.04948f9 3.14116f9 3.12580f9 … 8.60524f7 8.72117f7 8.30609f7
3.00483f9 2.98634f9 2.97949f9 8.54793f7 8.44453f7 8.84643f7
3.03344f9 3.00281f9 3.08042f9 8.42302f7 8.24331f7 8.63357f7
3.01833f9 3.02201f9 3.03027f9 8.60705f7 8.34739f7 8.34594f7
2.94578f9 2.99773f9 3.02747f9 8.04174f7 8.11727f7 8.16072f7
2.92023f9 2.92796f9 2.91467f9 … 7.80949f7 7.99472f7 7.81455f7
2.88401f9 2.95869f9 2.97999f9 7.91358f7 8.15653f7 7.97255f7
⋮ ⋱ ⋮
2.95995f9 2.90506f9 3.01727f9 … 8.05276f7 8.04053f7 8.1892f7
2.91456f9 2.92593f9 2.9851f9 7.91923f7 7.87905f7 7.93247f7
2.89651f9 2.81591f9 2.85835f9 7.95392f7 8.06558f7 7.94178f7
2.88585f9 2.86648f9 2.87025f9 8.20846f7 8.16997f7 8.20432f7
2.90665f9 2.81179f9 2.91064f9 8.04896f7 7.90742f7 8.22847f7
2.9013f9 2.93581f9 2.86071f9 … 8.0365f7 8.0742f7 8.21654f7
2.94684f9 2.95507f9 2.94315f9 7.8024f7 7.79604f7 7.95808f7Additionally, DimArrays have plot recipies that make visualization super easy. We can call Plots.jl's heatmap function to get a waterfall.
using Plots
fb.data |> heatmap
We can use this indexing to perform all sorts of analysis. Here we are using .. from IntervalSets.jl to create the range object.
fb.data[Freq = 1350..1450] |> heatmap
We can also index by time to get the nth integration
fb.data[Ti = 99] |> plot
We can also combine this indexing with Julia's builtin array operations such as looking at the time-averaged slice of spectrum. In this window, we get a nice look at the 21-centimeter line
using Statistics
dropdims(mean(fb.data,dims=Ti)[Freq = 1419..1422],dims=Ti) |> plot
Dedispersing
Let's look at some "realistic" data!
First we'll generate a fake dispersed pulse with a dispersion measure (DM) of 800 pc/cc over a frequency range of 1200-1500 MHz. We'll leave the other settings as is.
frb = pulse(800,1200,1500)A Filterbank file with 2048 time samples
----------
tsamp : 0.001000
Looking at the waterfall, we can see the frequency-integrated flux is spread out.
waterfall(frb)
We can apply the dedispersion transformation to get the real pulse shape
waterfall(dedisperse(frb,800))